% This file contains code to run simulations for the theory section of the 
% Handel, Kolstad, Spinniwijn paper on Adverse Selection and Information
% Frictions in Insurance Markets. 

% (Code below demonstrate two simulations, and additional simulations, 
%  assuming different parameters (levels of surplus, costs, and etc.)
%  have been omitted for brevity.)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%              First Simulation Set:                %%%
%%%                                                   %%%
%%% Similar to empirical environment                  %%%
%%% Low Surplus, High Frictions, Typical Costs        %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Distribution of expected total costs across population:
Costs = lognrnd(8,1.1,10000,1); 

% Generate SD of Costs for each level of overall costs
CostsSD = (4000 + 1.5*Costs);

% Distribution of Costs For Each Individual
Z = normrnd(0,1,1,50);
CostsTemp = (Costs*ones(1,50) + CostsSD*Z); 

% Final Distribution of Total Costs for Individuals
CostsDistribution = max(CostsTemp,0);

% Distribution of OOP costs in each of two plans
% Plan 1 is Free 
% Plan 2 has Deductible of 3000, Coinsurance 10%, OOPMax 7000
OOPCostsHDHP = zeros(10000,50); 
Temp1 = (CostsDistribution <=3000);
OOPCostsHDHP = OOPCostsHDHP + Temp1.*CostsDistribution;
Temp2 = (CostsDistribution >= 43000);
OOPCostsHDHP = OOPCostsHDHP + Temp2.*7000;
Temp3 = (1-Temp1).*(1-Temp2);
OOPCostsHDHP = OOPCostsHDHP + Temp3.*(3000 + 0.1*(CostsDistribution-3000));
 
% Compute Individual Expected OOP Costs in HDHP 
MeanOOPCostsHDHP = mean(OOPCostsHDHP,2);
 
% Compute Expected Firm Costs 
ExpectedCostsPPO = mean(CostsDistribution,2);
ExpectedCostsHDHP = mean(CostsDistribution - OOPCostsHDHP,2);

 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now, create distribution of frictions. In this first simulation, %
% these have high mean impact and substantial heterogeneity.       %
% Frictions here framed as mean impact on frictions of willingness %
% to pay for PPO.                                                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Frictions = normrnd(2500, 2000, 10000, 1); 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now, create Willingness to Pay and Surplus Values for Given Level of % 
% Risk Preferences. For this first simulation risk preferences         %
% are assumed to indicate low risk averaion, and low heterogeneity.    % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CARAPreferences = normrnd(0.0001, 0.00002, 10000,1);
CARAPreferences = CARAPreferences*ones(1,50);
 
% WTP For PPO For Perfectly Informed Consumer 
W = 100000; 
Utility = (W - OOPCostsHDHP);
UtilityHDHP = exp(-(CARAPreferences.*Utility));
CertaintyEquivalentHDHP = (-log(mean(UtilityHDHP,2))./mean(CARAPreferences,2));
 
WTPInformed = W - CertaintyEquivalentHDHP;    
 
% WTP For Consumers with Full Frictions 
WTPFullFrictions = WTPInformed + Frictions; 
WTPHalfFrictions = WTPInformed + 0.5*Frictions;
 
% HDHP Welfare Surplus (Above and Beyond Rel. Expected Costs Paid in HDHP)
Surplus = WTPInformed - (ExpectedCostsPPO - ExpectedCostsHDHP);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% RESULTS %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      %%% Simulation 1: Market with No Frictions %%%
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      
% Key curves for analysis: 
Quantity = (0:0.01:1);

DemandInformed = quantile(WTPInformed, 1-Quantity);
DemandFullFrictions = quantile(WTPFullFrictions, 1-Quantity);

% AC and MC curves conditional on demand curves
Quantiles = (.005:.01:.995)';

PPOMC = zeros(101,1);
HDHPMC = zeros(101,1);
NumbIn = zeros(101,1);

for n = 1:10000
    if WTPInformed(n,1) <= quantile(WTPInformed,.005)   
        PPOMC(101,1) = PPOMC(101,1) + ExpectedCostsPPO(n,1);
        HDHPMC(101,1) = HDHPMC(101,1) + ExpectedCostsHDHP(n,1);
        NumbIn(101,1) = NumbIn(101,1) + 1;
    end
    
    if WTPInformed(n,1) > quantile(WTPInformed,.995)  
        PPOMC(1,1) = PPOMC(1,1) + ExpectedCostsPPO(n,1);
        HDHPMC(1,1) = HDHPMC(1,1) + ExpectedCostsHDHP(n,1);
        NumbIn(1,1) = NumbIn(1,1) + 1;
    end
    for i = 1:99 
          if WTPInformed(n,1) <= quantile(WTPInformed,Quantiles(i+1,1)) && WTPInformed(n,1) > quantile(WTPInformed,Quantiles(i,1))  
                PPOMC(101-i,1) = PPOMC(101-i,1) + ExpectedCostsPPO(n,1);
                HDHPMC(101-i,1) = HDHPMC(101-i,1) + ExpectedCostsHDHP(n,1);
                NumbIn(101-i,1) = NumbIn(101-i,1) + 1;
          end
    end
end

PPOMCInformed = PPOMC./NumbIn;
HDHPMCInformed = HDHPMC./NumbIn; 

% Construct EFC Average Cost Curve
for i = 1:101
ACPPOEFC(i,1) = mean(PPOMCInformed(1:i));
end

% Construct HHW Average Cost Curve
for i = 1:101
ACPPOHHW(i,1) = mean(PPOMCInformed(1:i)) - mean(HDHPMCInformed(i:101));
end

% Compute Equilibrium Outcomes

    %%% EFC
    CrossingIndicator = 0;
    DemandSlot = 1:1:101;

    for i = 1:101
        if DemandInformed(1,i) <= ACPPOEFC(i,1)
           CrossingIndicator = CrossingIndicator + 1;
        end 
        if DemandInformed(1,i) > ACPPOEFC(i,1)
           CrossingIndicator = 0;
        end
        if CrossingIndicator == 1 
           QuantityEFC = (DemandSlot(1,i) - 1) / 100;
           PremiumEFC = ACPPOEFC(i,1);
        end 
    end

	%%% HHW 
    CrossingIndicator = 0;
    DemandSlot = 1:1:101;

    for i = 1:101
        if DemandInformed(1,i) <= ACPPOHHW(i,1)
           CrossingIndicator = CrossingIndicator + 1;
        end 
        if DemandInformed(1,i) > ACPPOHHW(i,1)
           CrossingIndicator = 0;
        end 
        if CrossingIndicator == 1 
           QuantityHHW = (DemandSlot(1,i) - 1) / 100;
           DeltaPremiumHHW = ACPPOHHW(i,1);
           PremiumHDHPHHW = mean(HDHPMCInformed(i:101));
        end 
    end

TotalSurplus = 0;
SurplusLostEFC = 0;
SurplusLostHHW = 0;

for i = 1:10000
    TotalSurplus = TotalSurplus + Surplus(i,1);
    
    if WTPInformed(i,1) < PremiumEFC 
        SurplusLostEFC = SurplusLostEFC + Surplus(i,1);
    end
    if WTPInformed(i,1) < DeltaPremiumHHW 
        SurplusLostHHW = SurplusLostHHW + Surplus(i,1);
    end
end 


      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      %%% Simulation 2: Market with Full Frictions %%%
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      
% Key curves for analysis: 
Quantity = (0:0.01:1);
DemandFullFrictions = quantile(WTPFullFrictions, 1-Quantity);

% AC and MC curves conditional on demand curves
Quantiles = (.005:.01:.995)';

PPOMC = zeros(101,1);
HDHPMC = zeros(101,1);
NumbIn = zeros(101,1);

for n = 1:10000
    if WTPFullFrictions(n,1) <= quantile(WTPFullFrictions,.005)   
        PPOMC(101,1) = PPOMC(101,1) + ExpectedCostsPPO(n,1);
        HDHPMC(101,1) = HDHPMC(101,1) + ExpectedCostsHDHP(n,1);
        NumbIn(101,1) = NumbIn(101,1) + 1;
    end
    if WTPFullFrictions(n,1) > quantile(WTPFullFrictions,.995)  
        PPOMC(1,1) = PPOMC(1,1) + ExpectedCostsPPO(n,1);
        HDHPMC(1,1) = HDHPMC(1,1) + ExpectedCostsHDHP(n,1);
        NumbIn(1,1) = NumbIn(1,1) + 1;
    end
    for i = 1:99 
          if WTPFullFrictions(n,1) <= quantile(WTPFullFrictions,Quantiles(i+1,1)) && WTPFullFrictions(n,1) > quantile(WTPFullFrictions,Quantiles(i,1))  
                PPOMC(101-i,1) = PPOMC(101-i,1) + ExpectedCostsPPO(n,1);
                HDHPMC(101-i,1) = HDHPMC(101-i,1) + ExpectedCostsHDHP(n,1);
                NumbIn(101-i,1) = NumbIn(101-i,1) + 1;
          end
    end
end

PPOMCFullFrictions = PPOMC./NumbIn;
HDHPMCFullFrictions = HDHPMC./NumbIn; 

% Construct EFC Average Cost Curve
for i = 1:101
ACPPOEFC(i,1) = mean(PPOMCFullFrictions(1:i));
end

% Construct HHW Average Cost Curve
for i = 1:101
ACPPOHHW(i,1) = mean(PPOMCFullFrictions(1:i)) - mean(HDHPMCFullFrictions(i:101));
end

% Compute Equilibrium Outcomes
    %%% EFC 
    CrossingIndicator = 0;
    DemandSlot = 1:1:101;

    for i = 1:101 
        if DemandFullFrictions(1,i) <= ACPPOEFC(i,1)
           CrossingIndicator = CrossingIndicator + 1;
        end 
        if DemandFullFrictions(1,i) > ACPPOEFC(i,1)
           CrossingIndicator = 0;
        end 
        if CrossingIndicator == 1 
           QuantityEFC = (DemandSlot(1,i) - 1) / 100;
           PremiumEFC = ACPPOEFC(i,1);
        end 
    end

    %%% HHW 
    CrossingIndicator = 0;
    DemandSlot = 1:1:101;

    for i = 1:101
        if DemandFullFrictions(1,i) <= ACPPOHHW(i,1)
           CrossingIndicator = CrossingIndicator + 1;
        end 
        if DemandFullFrictions(1,i) > ACPPOHHW(i,1)
           CrossingIndicator = 0;
        end 
        if CrossingIndicator == 1 
           QuantityHHW = (DemandSlot(1,i) - 1) / 100;
           DeltaPremiumHHW = ACPPOHHW(i,1);
           PremiumHDHPHHW = mean(HDHPMCFullFrictions(i:101));
        end 
    end

TotalSurplus =0;
SurplusLostEFC =0;
SurplusLostHHW = 0;

for i = 1:10000
    TotalSurplus = TotalSurplus + Surplus(i,1);
    if WTPFullFrictions(i,1) < PremiumEFC 
        SurplusLostEFC = SurplusLostEFC + Surplus(i,1);
    end
    if WTPFullFrictions(i,1) < DeltaPremiumHHW 
        SurplusLostHHW = SurplusLostHHW + Surplus(i,1);
    end
end 

      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%      
%%% Simulation 3: Market with Full Frictions and Half Risk Adjustment %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Beta = 0.5; 
MeanPPOCost = mean(ExpectedCostsPPO);
MeanHDHPCost = mean(ExpectedCostsHDHP);

% Key curves for analysis: 
Quantity = (0:0.01:1);
DemandFullFrictions = quantile(WTPFullFrictions, 1-Quantity);

% AC and MC curves conditional on demand curves

Quantiles = (.005:.01:.995)';

PPOMC = zeros(101,1);
HDHPMC = zeros(101,1);
NumbIn = zeros(101,1);

for n = 1:10000
    if WTPFullFrictions(n,1) <= quantile(WTPFullFrictions,.005)   
        PPOMC(101,1) = PPOMC(101,1) + (ExpectedCostsPPO(n,1)-Beta*(ExpectedCostsPPO(n,1)-MeanPPOCost));
        HDHPMC(101,1) = HDHPMC(101,1) + (ExpectedCostsHDHP(n,1)-Beta*(ExpectedCostsHDHP(n,1)-MeanHDHPCost));
        NumbIn(101,1) = NumbIn(101,1) + 1;
    end
    if WTPFullFrictions(n,1) > quantile(WTPFullFrictions,.995)  
        PPOMC(1,1) = PPOMC(1,1) + (ExpectedCostsPPO(n,1)-Beta*(ExpectedCostsPPO(n,1)-MeanPPOCost));
        HDHPMC(1,1) = HDHPMC(1,1) + (ExpectedCostsHDHP(n,1)-Beta*(ExpectedCostsHDHP(n,1)-MeanHDHPCost));
        NumbIn(1,1) = NumbIn(1,1) + 1;
    end
    for i = 1:99 
          if WTPFullFrictions(n,1) <= quantile(WTPFullFrictions,Quantiles(i+1,1)) && WTPFullFrictions(n,1) > quantile(WTPFullFrictions,Quantiles(i,1))  
                PPOMC(101-i,1) = PPOMC(101-i,1) + (ExpectedCostsPPO(n,1)--Beta*(ExpectedCostsPPO(n,1)-MeanPPOCost));
                HDHPMC(101-i,1) = HDHPMC(101-i,1) + (ExpectedCostsHDHP(n,1)-Beta*(ExpectedCostsHDHP(n,1)-MeanHDHPCost));
                NumbIn(101-i,1) = NumbIn(101-i,1) + 1;
          end
    end
end

PPOMCFullFrictions = PPOMC./NumbIn;
HDHPMCFullFrictions = HDHPMC./NumbIn; 

% Construct EFC Average Cost Curve
for i = 1:101
ACPPOEFC(i,1) = mean(PPOMCFullFrictions(1:i));
end

% Construct HHW Average Cost Curve
for i = 1:101
ACPPOHHW(i,1) = mean(PPOMCFullFrictions(1:i)) - mean(HDHPMCFullFrictions(i:101));
end

% Compute Equilibrium Outcomes
    %%% EFC 
    CrossingIndicator = 0;
    DemandSlot = 1:1:101;

    for i = 1:101
        if DemandFullFrictions(1,i) <= ACPPOEFC(i,1)
            CrossingIndicator = CrossingIndicator + 1;
        end 
        if DemandFullFrictions(1,i) > ACPPOEFC(i,1)
           CrossingIndicator = 0;
        end 
        if CrossingIndicator == 1 
           QuantityEFC = (DemandSlot(1,i) - 1) / 100;
           PremiumEFC = ACPPOEFC(i,1);
        end 
end

    %%% HHW 
    CrossingIndicator = 0;
    DemandSlot = 1:1:101;

    for i = 1:101
        if DemandFullFrictions(1,i) <= ACPPOHHW(i,1)
           CrossingIndicator = CrossingIndicator + 1;
        end 
        if DemandFullFrictions(1,i) > ACPPOHHW(i,1)
           CrossingIndicator = 0;
        end 
        if CrossingIndicator == 1 
           QuantityHHW = (DemandSlot(1,i) - 1) / 100;
           DeltaPremiumHHW = ACPPOHHW(i,1);
           PremiumHDHPHHW = mean(HDHPMCFullFrictions(i:101));
        end 
    end

TotalSurplus =0;
SurplusLostEFC =0;
SurplusLostHHW = 0;

for i = 1:10000
    TotalSurplus = TotalSurplus + Surplus(i,1);
    if WTPFullFrictions(i,1) < PremiumEFC 
        SurplusLostEFC = SurplusLostEFC + Surplus(i,1);
    end
    if WTPFullFrictions(i,1) < DeltaPremiumHHW 
        SurplusLostHHW = SurplusLostHHW + Surplus(i,1);
    end
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                 Second Simulation Set:               %%%
%%%                                                      %%%
%%% Med-High Surplus (mean and variance) High Frictions  %%%
%%% High Frictions towards PPO (mean and variance),      %%%
%%% Typical Costs (mean and variance)                    %%%
%%% No correlations between costs, frictions, risk pref  %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Distribution of expected total costs across population:
Costs = lognrnd(7.4,1.4,10000,1);
Costs = min(Costs,30000);
 
% Generate SD of Costs for each level of overall costs
CostsSD = (3000 + 1.2*Costs);

% Distribution of Costs For Each Individual
Z = normrnd(0,1,1,50);
CostsTemp = (Costs*ones(1,50) + CostsSD*Z); 

% Final Distribution of Total Costs for Individuals
CostsDistribution = max(CostsTemp,0);

% Distribution of OOP costs in each of two plans
% Plan 1 is Free 
% Plan 2 has Deductible of 3000, Coinsurance 10%, OOPMax 7000
OOPCostsHDHP = zeros(10000,50); 
Temp1 = (CostsDistribution <=3000);
OOPCostsHDHP = OOPCostsHDHP + Temp1.*CostsDistribution;
Temp2 = (CostsDistribution >= 43000);
OOPCostsHDHP = OOPCostsHDHP + Temp2.*7000;
Temp3 = (1-Temp1).*(1-Temp2);
OOPCostsHDHP = OOPCostsHDHP + Temp3.*(3000 + 0.1*(CostsDistribution-3000));
 
% Compute Individual Expected OOP Costs in HDHP 
MeanOOPCostsHDHP = mean(OOPCostsHDHP,2);
 
% Compute Expected Firm Costs 
ExpectedCostsPPO = mean(CostsDistribution,2);
ExpectedCostsHDHP = mean(CostsDistribution - OOPCostsHDHP,2);
 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now, create distribution of frictions. In this first simulation, these
% have high mean impact and substantial heterogeneity. 
% Frictions here framed as mean impact on frictions of willingness to pay
% for PPO. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Frictions = normrnd(2500, 2000,10000,1); 
Frictions1 = normrnd(2500,500,10000,1); 
Frictions2 = normrnd(0,2000,10000,1);
Frictions3 = normrnd(0,500,10000,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now, create Willingness to Pay and Surplus Values for Given Level of % 
% Risk Preferences. For this first simulation risk preferences         %
% are assumed to indicate low risk averaion, and low heterogeneity.    % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CARAPreferences = normrnd(0.001, 0.0003, 10000,1);
CARAPreferences = max(CARAPreferences,0.0001);
CARAPreferences = CARAPreferences*ones(1,50);
 
% WTP For PPO For Perfectly Informed Consumer 
W = 50000; 
Utility = (W - OOPCostsHDHP);
UtilityHDHP = exp(-(CARAPreferences.*Utility));
CertaintyEquivalentHDHP = (-log(mean(UtilityHDHP,2))./mean(CARAPreferences,2));
 
WTPInformed = W - CertaintyEquivalentHDHP;    
 
% WTP For Consumers with Full Frictions 
WTPFullFrictions = WTPInformed + Frictions;
WTPHalfFrictions = WTPInformed + 0.5*Frictions;
 
WTPFullFrictions1 = WTPInformed + Frictions1;
WTPHalfFrictions1 = WTPInformed + 0.5*Frictions1;
 
WTPFullFrictions2 = WTPInformed + Frictions2;
WTPHalfFrictions2 = WTPInformed + 0.5*Frictions2;
 
WTPFullFrictions3 = WTPInformed + Frictions3;
WTPHalfFrictions3 = WTPInformed + 0.5*Frictions3;
 
Surplus = WTPInformed - (ExpectedCostsPPO - ExpectedCostsHDHP);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% RESULTS %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Simulation 2-1: Market with Full Frictions (High Mean, High Var) %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Beta = 0; 
MeanPPOCost = mean(ExpectedCostsPPO);
MeanHDHPCost = mean(ExpectedCostsHDHP);

% Key curves for analysis: 
Quantity = (0:0.01:1);
DemandFullFrictions = quantile(WTPFullFrictions, 1-Quantity);

% AC and MC curves conditional on demand curves
Quantiles = (.005:.01:.995)';

PPOMC = zeros(101,1);
HDHPMC = zeros(101,1);
NumbIn = zeros(101,1);

for n = 1:10000
    if WTPFullFrictions(n,1) <= quantile(WTPFullFrictions,.005)   
        PPOMC(101,1) = PPOMC(101,1) + (ExpectedCostsPPO(n,1)-Beta*(ExpectedCostsPPO(n,1)-MeanPPOCost));
        HDHPMC(101,1) = HDHPMC(101,1) + (ExpectedCostsHDHP(n,1)-Beta*(ExpectedCostsHDHP(n,1)-MeanHDHPCost));
        NumbIn(101,1) = NumbIn(101,1) + 1;
    end
    if WTPFullFrictions(n,1) > quantile(WTPFullFrictions,.995)  
        PPOMC(1,1) = PPOMC(1,1) + (ExpectedCostsPPO(n,1)-Beta*(ExpectedCostsPPO(n,1)-MeanPPOCost));
        HDHPMC(1,1) = HDHPMC(1,1) + (ExpectedCostsHDHP(n,1)-Beta*(ExpectedCostsHDHP(n,1)-MeanHDHPCost));
        NumbIn(1,1) = NumbIn(1,1) + 1;
    end
    for i = 1:99 
          if WTPFullFrictions(n,1) <= quantile(WTPFullFrictions,Quantiles(i+1,1)) && WTPFullFrictions(n,1) > quantile(WTPFullFrictions,Quantiles(i,1))  
                PPOMC(101-i,1) = PPOMC(101-i,1) + (ExpectedCostsPPO(n,1)-Beta*(ExpectedCostsPPO(n,1)-MeanPPOCost));
                HDHPMC(101-i,1) = HDHPMC(101-i,1) + (ExpectedCostsHDHP(n,1)-Beta*(ExpectedCostsHDHP(n,1)-MeanHDHPCost));
                NumbIn(101-i,1) = NumbIn(101-i,1) + 1;
          end
    end
end

PPOMCFullFrictions = PPOMC./NumbIn;
HDHPMCFullFrictions = HDHPMC./NumbIn; 

% Construct EFC Average Cost Curve
for i = 1:101
ACPPOEFC(i,1) = mean(PPOMCFullFrictions(1:i)) - mean(HDHPMCFullFrictions(1:i));
end

% Construct HHW Average Cost Curve
for i = 1:101
ACPPOHHW(i,1) = mean(PPOMCFullFrictions(1:i)) - mean(HDHPMCFullFrictions(i:101));
end

% Compute Equilibrium Outcomes
    %%% EFC
    CrossingIndicator = 0;
    DemandSlot = 1:1:101;

    for i = 1:101
        if DemandFullFrictions(1,i) <= ACPPOEFC(i,1)
           CrossingIndicator = CrossingIndicator + 1;
        end 
        if DemandFullFrictions(1,i) > ACPPOEFC(i,1)
           CrossingIndicator = 0;
        end 
        if CrossingIndicator == 1 
           QuantityEFC = (DemandSlot(1,i) - 1) / 100;
           PremiumEFC = ACPPOEFC(i,1);
        end 
    end

    %%% HHW 
    CrossingIndicator = 0;
    DemandSlot = 1:1:101;

    for i = 1:101
        if DemandFullFrictions(1,i) <= ACPPOHHW(i,1)
           CrossingIndicator = CrossingIndicator + 1;
        end 
        if DemandFullFrictions(1,i) > ACPPOHHW(i,1)
           CrossingIndicator = 0;
        end 
        if CrossingIndicator == 1 
           QuantityHHW = (DemandSlot(1,i) - 1) / 100;
           DeltaPremiumHHW = ACPPOHHW(i,1);
           PremiumHDHPHHW = mean(HDHPMCFullFrictions(i:101));
        end 
    end

TotalSurplus = 0;
SurplusLostEFC = 0;
SurplusLostHHW = 0;

for i = 1:10000
    TotalSurplus = TotalSurplus + Surplus(i,1);
    
    if WTPFullFrictions(i,1) < PremiumEFC 
        SurplusLostEFC = SurplusLostEFC + Surplus(i,1);
    end
    if WTPFullFrictions(i,1) < DeltaPremiumHHW 
        SurplusLostHHW = SurplusLostHHW + Surplus(i,1);
    end
end 

% Compute Proportion of Consumers with Frictions Making Choices
% Where they make the wrong choice because of frictions given the
% equilibrium premium. 
MistakesLowEFC = 0;
MistakesHighEFC = 0;

for i = 1:10000
    if WTPInformed(i,1) <= PremiumEFC && PremiumEFC <= WTPFullFrictions(i,1)
       MistakesHighEFC = MistakesHighEFC + 1;    
    end
    if WTPFullFrictions(i,1) <= PremiumEFC && PremiumEFC <= WTPInformed(i,1)
       MistakesLowEFC = MistakesLowEFC + 1;    
    end
end

MistakesLow = 0;
MistakesHigh = 0;

for i = 1:10000
    if WTPInformed(i,1) <= DeltaPremiumHHW && DeltaPremiumHHW <= WTPFullFrictions(i,1)
       MistakesHigh = MistakesHigh + 1;    
    end
    if WTPFullFrictions(i,1) <= DeltaPremiumHHW && DeltaPremiumHHW <= WTPInformed(i,1)
       MistakesLow = MistakesLow + 1;    
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Simulation 2-2: Market with Half Frictions (High Mean, High Var)  %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Beta = 0; 
MeanPPOCost = mean(ExpectedCostsPPO);
MeanHDHPCost = mean(ExpectedCostsHDHP);

% Key curves for analysis: 
Quantity = (0:0.01:1);
DemandHalfFrictions = quantile(WTPHalfFrictions, 1-Quantity);

% AC and MC curves conditional on demand curves
Quantiles = (.005:.01:.995)';

PPOMC = zeros(101,1);
HDHPMC = zeros(101,1);
NumbIn = zeros(101,1);

for n = 1:10000
    if WTPHalfFrictions(n,1) <= quantile(WTPHalfFrictions,.005)   
        PPOMC(101,1) = PPOMC(101,1) + (ExpectedCostsPPO(n,1)-Beta*(ExpectedCostsPPO(n,1)-MeanPPOCost));
        HDHPMC(101,1) = HDHPMC(101,1) + (ExpectedCostsHDHP(n,1)-Beta*(ExpectedCostsHDHP(n,1)-MeanHDHPCost));
        NumbIn(101,1) = NumbIn(101,1) + 1;
    end
    
    if WTPHalfFrictions(n,1) > quantile(WTPHalfFrictions,.995)  
        PPOMC(1,1) = PPOMC(1,1) + (ExpectedCostsPPO(n,1)-Beta*(ExpectedCostsPPO(n,1)-MeanPPOCost));
        HDHPMC(1,1) = HDHPMC(1,1) + (ExpectedCostsHDHP(n,1)-Beta*(ExpectedCostsHDHP(n,1)-MeanHDHPCost));
        NumbIn(1,1) = NumbIn(1,1) + 1;
    end
    
    for i = 1:99 
          if WTPHalfFrictions(n,1) <= quantile(WTPHalfFrictions,Quantiles(i+1,1)) && WTPHalfFrictions(n,1) > quantile(WTPHalfFrictions,Quantiles(i,1))  
                PPOMC(101-i,1) = PPOMC(101-i,1) + (ExpectedCostsPPO(n,1)-Beta*(ExpectedCostsPPO(n,1)-MeanPPOCost));
                HDHPMC(101-i,1) = HDHPMC(101-i,1) + (ExpectedCostsHDHP(n,1)-Beta*(ExpectedCostsHDHP(n,1)-MeanHDHPCost));
                NumbIn(101-i,1) = NumbIn(101-i,1) + 1;
          end
    end
end

PPOMCHalfFrictions = PPOMC./NumbIn;
HDHPMCHalfFrictions = HDHPMC./NumbIn; 

% Construct EFC Average Cost Curve
for i = 1:101
ACPPOEFC(i,1) = mean(PPOMCHalfFrictions(1:i)) - mean(HDHPMCHalfFrictions(1:i));
end

% Construct HHW Average Cost Curve
for i = 1:101
ACPPOHHW(i,1) = mean(PPOMCHalfFrictions(1:i)) - mean(HDHPMCHalfFrictions(i:101));
end

% Compute Equilibrium Outcomes
    %%% EFC
    CrossingIndicator = 0;
    DemandSlot = 1:1:101;

    for i = 1:101
        if DemandHalfFrictions(1,i) <= ACPPOEFC(i,1)
           CrossingIndicator = CrossingIndicator + 1;
        end 
        if DemandHalfFrictions(1,i) > ACPPOEFC(i,1)
           CrossingIndicator = 0;
        end 
        if CrossingIndicator == 1 
           QuantityEFC = (DemandSlot(1,i) - 1) / 100;
           PremiumEFC = ACPPOEFC(i,1);
        end 
    end

    %%% HHW 
    CrossingIndicator = 0;
    DemandSlot = 1:1:101;

    for i = 1:101
        if DemandHalfFrictions(1,i) <= ACPPOHHW(i,1)
           CrossingIndicator = CrossingIndicator + 1;
        end 
        if DemandHalfFrictions(1,i) > ACPPOHHW(i,1)
           CrossingIndicator = 0;
        end 
        if CrossingIndicator == 1 
           QuantityHHW = (DemandSlot(1,i) - 1) / 100;
           DeltaPremiumHHW = ACPPOHHW(i,1);
           PremiumHDHPHHW = mean(HDHPMCHalfFrictions(i:101));
        end 
    end

TotalSurplus =0;
SurplusLostEFC =0;
SurplusLostHHW = 0;

for i = 1:10000
    TotalSurplus = TotalSurplus + Surplus(i,1);
    if WTPHalfFrictions(i,1) < PremiumEFC 
        SurplusLostEFC = SurplusLostEFC + Surplus(i,1);
    end
    if WTPHalfFrictions(i,1) < DeltaPremiumHHW 
        SurplusLostHHW = SurplusLostHHW + Surplus(i,1);
    end
end 

% Compute Proportion of Consumers with Frictions Making Choices
% Where they make the wrong choice because of frictions given the
% equilibrium premium. 

MistakesLowEFC = 0;
MistakesHighEFC = 0;

for i = 1:10000
    if WTPInformed(i,1) <= PremiumEFC && PremiumEFC <= WTPHalfFrictions(i,1)
       MistakesHighEFC = MistakesHighEFC + 1;    
    end
    if WTPHalfFrictions(i,1) <= PremiumEFC && PremiumEFC <= WTPInformed(i,1)
       MistakesLowEFC = MistakesLowEFC + 1;    
    end
end

MistakesLow = 0;
MistakesHigh = 0;

for i = 1:10000
    if WTPInformed(i,1) <= DeltaPremiumHHW && DeltaPremiumHHW <= WTPHalfFrictions(i,1)
       MistakesHigh = MistakesHigh + 1;    
    end
    if WTPHalfFrictions(i,1) <= DeltaPremiumHHW && DeltaPremiumHHW <= WTPInformed(i,1)
       MistakesLow = MistakesLow + 1;
    end
end


        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%% Simulation 2-3: Market with No Frictions %%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Beta = 0; 
MeanPPOCost = mean(ExpectedCostsPPO);
MeanHDHPCost = mean(ExpectedCostsHDHP);

% Key curves for analysis: 
Quantity = (0:0.01:1);
DemandInformed = quantile(WTPInformed, 1-Quantity);

% Now create AC and MC curves conditional on demand curves
Quantiles = (.005:.01:.995)';

PPOMC = zeros(101,1);
HDHPMC = zeros(101,1);
NumbIn = zeros(101,1);

for n = 1:10000
    if WTPInformed(n,1) <= quantile(WTPInformed,.005)   
        PPOMC(101,1) = PPOMC(101,1) + (ExpectedCostsPPO(n,1)-Beta*(ExpectedCostsPPO(n,1)-MeanPPOCost));
        HDHPMC(101,1) = HDHPMC(101,1) + (ExpectedCostsHDHP(n,1)-Beta*(ExpectedCostsHDHP(n,1)-MeanHDHPCost));
        NumbIn(101,1) = NumbIn(101,1) + 1;
    end
    
    if WTPInformed(n,1) > quantile(WTPInformed,.995)  
        PPOMC(1,1) = PPOMC(1,1) + (ExpectedCostsPPO(n,1)-Beta*(ExpectedCostsPPO(n,1)-MeanPPOCost));
        HDHPMC(1,1) = HDHPMC(1,1) + (ExpectedCostsHDHP(n,1)-Beta*(ExpectedCostsHDHP(n,1)-MeanHDHPCost));
        NumbIn(1,1) = NumbIn(1,1) + 1;
    end
    for i = 1:99 
          if WTPInformed(n,1) <= quantile(WTPInformed,Quantiles(i+1,1)) && WTPInformed(n,1) > quantile(WTPInformed,Quantiles(i,1))  
                PPOMC(101-i,1) = PPOMC(101-i,1) + (ExpectedCostsPPO(n,1)-Beta*(ExpectedCostsPPO(n,1)-MeanPPOCost));
                HDHPMC(101-i,1) = HDHPMC(101-i,1) + (ExpectedCostsHDHP(n,1)-Beta*(ExpectedCostsHDHP(n,1)-MeanHDHPCost));
                NumbIn(101-i,1) = NumbIn(101-i,1) + 1;
          end
    end
end

PPOMCInformed = PPOMC./NumbIn;
HDHPMCInformed = HDHPMC./NumbIn; 

% Construct EFC Average Cost Curve
for i = 1:101
ACPPOEFC(i,1) = mean(PPOMCInformed(1:i)) - mean(HDHPMCInformed(1:i));
end

% Construct HHW Average Cost Curve
for i = 1:101
ACPPOHHW(i,1) = mean(PPOMCInformed(1:i)) - mean(HDHPMCInformed(i:101));
end

% Compute Equilibrium Outcomes
    %%% EFC
    CrossingIndicator = 0;
    DemandSlot = 1:1:101;

    for i = 1:101
        if DemandInformed(1,i) <= ACPPOEFC(i,1)
           CrossingIndicator = CrossingIndicator + 1;
        end 
        if DemandInformed(1,i) > ACPPOEFC(i,1)
           CrossingIndicator = 0;
        end 
        if CrossingIndicator == 1 
           QuantityEFC = (DemandSlot(1,i) - 1) / 100;
           PremiumEFC = ACPPOEFC(i,1);
        end 
    end

    %%% EFC
    CrossingIndicator = 0;
    DemandSlot = 1:1:101;

    for i = 1:101
        if DemandInformed(1,i) <= ACPPOHHW(i,1)
           CrossingIndicator = CrossingIndicator + 1;
        end 
        if DemandInformed(1,i) > ACPPOHHW(i,1)
           CrossingIndicator = 0;
        end 
        if CrossingIndicator == 1 
           QuantityHHW = (DemandSlot(1,i) - 1) / 100;
           DeltaPremiumHHW = ACPPOHHW(i,1);
           PremiumHDHPHHW = mean(HDHPMCInformed(i:101));
        end 
    end

TotalSurplus =0;
SurplusLostEFC =0;
SurplusLostHHW = 0;

for i = 1:10000
    TotalSurplus = TotalSurplus + Surplus(i,1);
    if WTPInformed(i,1) < PremiumEFC 
        SurplusLostEFC = SurplusLostEFC + Surplus(i,1);
    end
    if WTPInformed(i,1) < DeltaPremiumHHW 
        SurplusLostHHW = SurplusLostHHW + Surplus(i,1);
    end
end 

% Compute Proportion of Consumers with Frictions Making Choices
% Where they make the wrong choice because of frictions given the
% equilibrium premium. 

MistakesLowEFC = 0;
MistakesHighEFC = 0;

for i = 1:10000 
    if WTPInformed(i,1) <= PremiumEFC && PremiumEFC <= WTPInformed(i,1)
       MistakesHighEFC = MistakesHighEFC + 1;    
    end
    if WTPInformed(i,1) <= PremiumEFC && PremiumEFC <= WTPInformed(i,1)
       MistakesLowEFC = MistakesLowEFC + 1;     
    end
end

MistakesLow = 0;
MistakesHigh = 0;

for i = 1:10000 
    if WTPInformed(i,1) <= DeltaPremiumHHW && DeltaPremiumHHW <= WTPInformed(i,1)
       MistakesHigh = MistakesHigh + 1;    
    end
    if WTPInformed(i,1) <= DeltaPremiumHHW && DeltaPremiumHHW <= WTPInformed(i,1)
       MistakesLow = MistakesLow + 1;    
    end
end